﻿using csharpcl;
using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.IO;
using System.Linq;
using System.Threading.Tasks;

namespace CL.Demo
{
    public class Program
    {
        public const double Double_Zero = 1E-8D;
        public const float Float_Zero = 1E-8F;

        #region 矩阵函数

        private static void WriteFile(double[] cx, int Ndim, int Mdim, string file, string headertext)
        {
            if (File.Exists(file))
                File.Delete(file);
            using (var w = new StreamWriter(File.OpenWrite(file)))
            {
                w.WriteLine(headertext);
                w.WriteLine();
                w.WriteLine("矩阵大小:{0}x{1}", Ndim, Mdim);
                w.WriteLine();
                for (int n = 0; n < Ndim; n++)
                {
                    for (int m = 0; m < Mdim; m++)
                    {
                        int index = n * Mdim + m;
                        w.Write(string.Format("{0}\t", cx[index].ToString("F6")));
                    }
                    w.WriteLine();
                }
                w.Flush();
            }
        }

        private static void WriteFileTolorance(Tuple<int, double, double, double>[] tpls, string file, string headertext)
        {
            if (File.Exists(file))
                File.Delete(file);
            using (var w = new StreamWriter(File.OpenWrite(file)))
            {
                w.WriteLine(headertext);
                w.WriteLine();
                w.WriteLine("索引\t误差\t\t值1\t\t值2");
                w.WriteLine("----\t---\t\t---\t\t---");
                for (int j = 0; j < tpls.Length; j++)
                {
                    w.WriteLine(string.Format("{0}\t{1:F8}\t{2:F8}\t{3:F8}", tpls[j].Item1, tpls[j].Item2, tpls[j].Item3, tpls[j].Item4));
                }
                w.Flush();
            }
        }

        private static int GetIndex(int i, int j, int Ndim, int Mdim)
        {
            return i * Mdim + j;
        }

        private static void PrintMatrix(float[] matrix, int Ndim, int Mdim)
        {
            if (matrix == null)
            {
                Console.WriteLine("NULL Matrix");
                return;
            }
            for (int j = 0; j < Mdim; j++)
            {
                for (int i = 0; i < Ndim; i++)
                    Console.Write("{0}\t", matrix[j * Ndim + i].ToString("F6"));
                Console.WriteLine();
            }
        }

        private static void PrintMatirxNoZero(float[] matrix, int Ndim, int Mdim, float zero = Float_Zero)
        {
            for (int i = 0; i < Ndim; i++)
            {
                for (int j = 0; j < Mdim; j++)
                {
                    if (Math.Abs(matrix[i * Mdim + j]) > zero)
                        Console.Write("{0}\t", matrix[i * Mdim + j]);
                }
                //Console.WriteLine();
            }
        }

        private static Tuple<int, float, float, float>[] BiggerZero(float[] A, float[] B, float zero = Float_Zero)
        {
            var lis = new List<Tuple<int, float, float, float>>();
            for (int i = 0; i < A.Length; i++)
            {
                float d = A[i] - B[i];
                if (Math.Abs(d) > zero)
                {
                    lis.Add(Tuple.Create(i, d, A[i], B[i]));
                }
            }
            return lis.ToArray();
        }

        private static Tuple<int, double, double, double>[] BiggerZero(double[] A, double[] B, double zero = Double_Zero)
        {
            var lis = new List<Tuple<int, double, double, double>>();
            for (int i = 0; i < A.Length; i++)
            {
                double d = A[i] - B[i];
                if (Math.Abs(d) > zero)
                {
                    lis.Add(Tuple.Create(i, d, A[i], B[i]));
                }
            }
            return lis.ToArray();
        }

        /// <summary>
        /// 矩阵相减
        /// </summary>
        /// <param name="A"></param>
        /// <param name="B"></param>
        /// <returns></returns>
        private static float[] SubMatrix(float[] A, float[] B)
        {
            var C = new float[A.Length];
            for (int i = 0; i < A.Length; i++)
            {
                C[i] = A[i] - B[i];
            }
            return C;
        }

        private static double[] SubMatrix(double[] A, float[] B)
        {
            var C = new double[A.Length];
            for (int i = 0; i < A.Length; i++)
            {
                C[i] = A[i] - B[i];
            }
            return C;
        }

        /// <summary>
        /// 构建单位矩阵
        /// </summary>
        /// <param name="A"></param>
        /// <param name="Ndim"></param>
        /// <param name="Mdim"></param>
        private static void UnitMatrix(float[] A, int Ndim, int Mdim)
        {
            for (int i = 0; i < Ndim; i++)
            {
                for (int j = 0; j < Mdim; j++)
                {
                    int ii = i * Mdim;
                    if (i == j)
                        A[ii + j] = 1;
                    else
                        A[ii + j] = 0;
                }
            }
        }

        private static void UnitMatrix(double[] A, int Ndim, int Mdim)
        {
            for (int i = 0; i < Ndim; i++)
            {
                for (int j = 0; j < Mdim; j++)
                {
                    int ii = i * Mdim;
                    if (i == j)
                        A[ii + j] = 1;
                    else
                        A[ii + j] = 0;
                }
            }
        }

        private static void NumberMatrix(float[] A, int Ndim = 0, int Mdim = 0)
        {
            for (int i = 0; i < A.Length; i++)
            {
                A[i] = 1 + i * 0.00001f;
            }
        }

        private static void NumberMatrix(double[] A, int Ndim = 0, int Mdim = 0)
        {
            for (int i = 0; i < A.Length; i++)
            {
                A[i] = 1 + 0.00001d * i;
            }
        }

        private static void NumberMatrix2(double[] A, int Ndim = 0, int Mdim = 0)
        {
            for (int i = 0; i < A.Length; i++)
            {
                A[i] = i + 0.0000001d * i;
            }
        }

        private static void SpecMatrix(float[] A, int Ndim = 0, int Mdim = 0)
        {
            for (int j = 0; j < Mdim; j++)
            {
                for (int i = 0; i < Ndim; i++)
                {
                    int ii = j * Ndim;
                    A[ii + i] = i + 1;
                }
            }
        }

        #endregion 矩阵函数

        #region 矩阵乘法

        public const string CL_MulMatrix_Float = @"
__kernel void matrix_mult(
	const int Ndim,
	const int Mdim,
	const int Pdim,
	__global const float* A,
	__global const float* B,
	__global float* C)
{
	int i = get_global_id(0);
	int j = get_global_id(1);

	int k;
	float tmp;

	if ((i < Ndim) && (j < Mdim)) {
		tmp = 0.0;
		for (k = 0; k < Pdim; k++)
			tmp += A[i * Pdim + k] * B[k * Mdim + j];
		C[i * Mdim + j] = tmp;
	}
}
";

        public const string CL_MulMatrix_Double = @"
__kernel void matrix_mult(
	const int Ndim,
	const int Mdim,
	const int Pdim,
	__global const double* A,
	__global const double* B,
	__global double* C)
{
	int i = get_global_id(0);
	int j = get_global_id(1);

	int k;
	double tmp;

	if ((i < Ndim) && (j < Mdim)) {
		tmp = 0.0;
		for (k = 0; k < Pdim; k++)
			tmp += A[i * Pdim + k] * B[k * Mdim + j];
		C[i * Mdim + j] = tmp;
	}
}
";

        public static void TestMulMatrix()
        {
            Console.Title = "矩阵计算测试";
            Console.WriteLine("矩阵计算测试");

            int x = 11;
            Introp.Add(77, 11, ref x);

            string file = "./addVec.cl";
            Console.WriteLine("输入要计算的方阵大小:...");
            int num = 1000;
            while (true)
            {
                var str = Console.ReadLine();
                if (int.TryParse(str, out num))
                {
                    if (num <= 0)
                    {
                        Console.WriteLine("无效的方阵大小");
                    }
                    else
                        break;
                }
                Console.WriteLine("输入计算的方阵大小:...");
            }

            int Ndim = num;
            int Mdim = num;
            int Pdim = num;
            int szA = Ndim * Pdim;
            int szB = Pdim * Mdim;
            int szC = Ndim * Mdim;

            Console.WriteLine("----------计算的矩阵大小----------");
            Console.WriteLine("A:{0}x{1},B:{1}x{2},C:{0}x{2}", Ndim, Pdim, Mdim);
            Console.WriteLine("A x B = C");
            Console.WriteLine("---------------------------------");
            Console.WriteLine();

            double[] A = new double[szA];
            double[] B = new double[szB];
            double[] C, cx, cx2;

            //初始化矩阵
            NumberMatrix2(A, Ndim, Pdim);
            NumberMatrix(B, Pdim, Mdim);
            //UnitMatrix(B, Pdim, Mdim);
            //NumberMatrix(B, Pdim, Mdim);
            //ZeroMatrix(C);
            //ZeroMatrix(_outCX);
            //ZeroMatrix(_outCX2);

            double[] ts = new double[10];
            var st = new Stopwatch();
            Console.WriteLine("----------GPU计算开始----------");
            st.Start();
            //执行计算
            C = MulByGPU(A, B, Ndim, Pdim, Mdim);

            st.Stop();
            Console.WriteLine("----------GPU计算结束----------");
            ts[0] = st.ElapsedMilliseconds;
            Console.WriteLine("GPU耗时(ms):\t" + ts[0]);
            Console.WriteLine();
            Console.WriteLine("----------CPU1并行计算开始----------");
            st.Restart();

            #region CPU并行计算

            cx = MulByCPUAsync(A, B, Ndim, Pdim, Mdim);

            #endregion CPU并行计算

            st.Stop();
            ts[1] = st.ElapsedMilliseconds;
            Console.WriteLine("----------CPU1并行计算结束----------");
            Console.WriteLine();
            Console.WriteLine("----------CPU2并行计算开始----------");
            st.Restart();

            #region CPU并行计算

            cx2 = MulByCPUAsync(A, B, Ndim, Pdim, Mdim);

            #endregion CPU并行计算

            st.Stop();
            ts[2] = st.ElapsedMilliseconds;
            Console.WriteLine("----------CPU2并行计算结束----------");

#if CPU_SYNC
            Console.WriteLine("----------CPU串行计算开始----------");
            st.Restart();

            #region CPU串行计算

            Mul(A, B, _outCX2, Ndim, Pdim, Mdim);

            #endregion CPU串行计算

            st.Stop();
            t2 = st.ElapsedMilliseconds;
            Console.WriteLine("----------CPU串行计算结束----------");
#endif
            Console.WriteLine();
            Console.WriteLine("************************");
            Console.WriteLine("GPU      耗时(ms):\t" + ts[0]);
            Console.WriteLine("CPU 1并行耗时(ms):\t" + ts[1]);
            Console.WriteLine("CPU 2并行耗时(ms):\t" + ts[2]);
            Console.WriteLine("************************");

            Console.WriteLine();

            Tuple<int, double, double, double>[] tpls = BiggerZero(cx, C);
            Tuple<int, double, double, double>[] tpls2 = BiggerZero(cx, cx2);

            Console.WriteLine("误差大于【{0:E}】", Double_Zero);
            Console.WriteLine("[CPU并行比较GPU]        数量：" + tpls.Count());
            Console.WriteLine("[CPU1并行比较CPU2并行]  数量：" + tpls2.Count());

            Console.WriteLine();

            Console.WriteLine("将【CPU对比GPU】结果写入文件");
            WriteFileTolorance(tpls, "00-对比结果-CPU-GPU.txt", string.Format("CPU对比GPU,误差大于:{0:E}", Double_Zero));

            Console.WriteLine("将【CPU1对比CPU2】结果写入文件");
            WriteFileTolorance(tpls2, "00-对比结果-CPU1-CPU2.txt", string.Format("CPU1对比CPU2,误差大于:{0:E}", Double_Zero));

            Console.WriteLine("将【输入矩阵A】结果写入文件");
            WriteFile(A, Ndim, Pdim, "00-输入矩阵-A.txt", "输入矩阵A");

            Console.WriteLine("将【输入矩阵B】结果写入文件");
            WriteFile(B, Pdim, Mdim, "00-输入矩阵-B.txt", "输入矩阵B");

            Console.WriteLine("将【GPU计算结果】结果写入文件");
            WriteFile(C, Ndim, Mdim, "00-计算结果-GPU-01.txt", "GPU 1计算结果，耗时(ms):" + ts[0].ToString());

            Console.WriteLine("将【CPU1计算结果】结果写入文件");
            WriteFile(cx, Ndim, Mdim, "00-计算结果-CPU-01.txt", "CPU 1并行计算结果，耗时(ms):" + ts[1].ToString());

            Console.WriteLine("将【CPU2计算结果】结果写入文件");
            WriteFile(cx2, Ndim, Mdim, "00-计算结果-CPU-02.txt", "CPU 2并行计算结果，耗时(ms):" + ts[2].ToString());
        }

        #region CPU矩阵乘法

        /// <summary>
        /// CPU串行的矩阵乘法
        /// </summary>
        /// <param name="A"></param>
        /// <param name="B"></param>
        /// <param name="Ndim"></param>
        /// <param name="Pdim"></param>
        /// <param name="Mdim"></param>
        /// <returns></returns>
        private static float[] MulByCPUSync(float[] A, float[] B, int Ndim, int Pdim, int Mdim)
        {
            var C = new float[Ndim * Mdim];
            //计算矩阵乘法
            for (int i = 0; i < Ndim; i++)
            {
                for (int j = 0; j < Mdim; j++)
                {
                    float sum = 0;
                    for (int k = 0; k < Pdim; k++)
                    {
                        sum += A[i * Pdim + k] * B[k * Pdim + j];
                    }
                    C[i * Mdim + j] = sum;
                }
            }
            return C;
        }

        private static double[] MulByCPUSync(double[] A, double[] B, int Ndim, int Pdim, int Mdim)
        {
            var C = new double[Ndim * Mdim];
            //计算矩阵乘法
            for (int i = 0; i < Ndim; i++)
            {
                for (int j = 0; j < Mdim; j++)
                {
                    double sum = 0;
                    for (int k = 0; k < Pdim; k++)
                    {
                        sum += A[i * Pdim + k] * B[k * Pdim + j];
                    }
                    C[i * Mdim + j] = sum;
                }
            }
            return C;
        }

        /// <summary>
        /// CPU按行并行的矩阵乘法
        /// </summary>
        /// <param name="A"></param>
        /// <param name="B"></param>
        /// <param name="Ndim"></param>
        /// <param name="Pdim"></param>
        /// <param name="Mdim"></param>
        /// <returns></returns>
        private static float[] MulByCPUAsync(float[] A, float[] B, int Ndim, int Pdim, int Mdim)
        {
            var C = new float[Ndim * Mdim];
            //计算矩阵乘法
            //for (int i = 0; i < Ndim; i++)
            Parallel.For(0, Ndim, i =>
            {
                for (int j = 0; j < Mdim; j++)
                {
                    float sum = 0;
                    for (int k = 0; k < Pdim; k++)
                    {
                        sum += A[i * Pdim + k] * B[k * Pdim + j];
                    }
                    C[i * Mdim + j] = sum;
                }
            });
            return C;
        }

        private static double[] MulByCPUAsync(double[] A, double[] B, int Ndim, int Pdim, int Mdim)
        {
            var C = new double[Ndim * Mdim];
            //计算矩阵乘法
            //for (int i = 0; i < Ndim; i++)
            Parallel.For(0, Ndim, i =>
            {
                for (int j = 0; j < Mdim; j++)
                {
                    double sum = 0;
                    for (int k = 0; k < Pdim; k++)
                    {
                        sum += A[i * Pdim + k] * B[k * Pdim + j];
                    }
                    C[i * Mdim + j] = sum;
                }
            });
            return C;
        }

        #endregion CPU矩阵乘法

        #region GPU矩阵乘法

        /// <summary>
        /// GPU计算矩阵乘法
        /// </summary>
        /// <param name="A"></param>
        /// <param name="B"></param>
        /// <param name="Ndim"></param>
        /// <param name="Pdim"></param>
        /// <param name="Mdim"></param>
        /// <returns></returns>
        private static float[] MulByGPU(float[] A, float[] B, int Ndim, int Pdim, int Mdim)
        {
            float[] C = new float[Ndim * Mdim];

            ulong pid = 0;
            ulong did = 0;
            IntPtr ctx = IntPtr.Zero;
            IntPtr que = IntPtr.Zero;
            IntPtr prg = IntPtr.Zero;
            IntPtr kel = IntPtr.Zero;

            var pids = Introp.GetPlatformIds();
            //if (pids == null || pids.Length == 0)
            //{
            //    Console.WriteLine("zero platform");
            //    Console.ReadLine();
            //    return;
            //}
            var dic = new Dictionary<ulong, ulong[]>();
            for (int i = 0; i < pids.Length; i++)
            {
                dic.Add(pids[i], Introp.GetDeviceIds(pids[i]));
            }
            //foreach (var item in dic)
            //{
            //    Console.WriteLine(Introp.GetPlatformName(item.Key));
            //    if (item.Value == null || item.Value.Length == 0)
            //        continue;
            //    for (int i = 0; i < item.Value.Length; i++)
            //    {
            //        Console.WriteLine("\t" + Introp.GetDeviceName(item.Value[i]));
            //    }
            //}
            pid = dic.Keys.ToArray()[0];
            did = dic[pid][0];
            Console.WriteLine("选择的平台：" + Introp.GetPlatformName(pid));
            Console.WriteLine("选择的设备：" + Introp.GetDeviceName(did));
            ctx = Introp.CreateContext(did);
            que = Introp.CreateCommandQueue(did, ctx);
            string err;
            //使用cl文件
            //prg = Introp.CreateProgramFile(did, ctx, file, out err);
            //使用内置字符串
            prg = Introp.CreateProgramCode(did, ctx, CL_MulMatrix_Float, out err);
            kel = Introp.CreateKernal(prg, "matrix_mult");

            IntPtr pA = Introp.CreateClInputBuffer(ctx, A);
            IntPtr pB = Introp.CreateClInputBuffer(ctx, B);
            IntPtr pC = Introp.CreateClOutputBuffer(ctx, C);

            Introp.SetKernalParam(kel, 0, Ndim);
            Introp.SetKernalParam(kel, 1, Mdim);
            Introp.SetKernalParam(kel, 2, Pdim);
            Introp.SetKernalParam(kel, 3, pA);
            Introp.SetKernalParam(kel, 4, pB);
            Introp.SetKernalParam(kel, 5, pC);
            Int64[] glbal = new Int64[]
            {
                Ndim,Mdim
            };

            IntPtr pevent = Introp.Exec_CommandQueue(que, kel, 2, glbal);

            Introp.Exec_CommandQueueReadBuffer(que, pC, C);

            Introp.ReleaseMem(pA);
            Introp.ReleaseMem(pB);
            Introp.ReleaseMem(pC);
            Introp.ReleaseCommandQueue(que);
            Introp.ReleaseProgram(prg);
            Introp.ReleaseContext(ctx);

            return C;
        }

        private static double[] MulByGPU(double[] A, double[] B, int Ndim, int Pdim, int Mdim)
        {
            double[] C = new double[Ndim * Mdim];

            ulong pid = 0;
            ulong did = 0;
            IntPtr ctx = IntPtr.Zero;
            IntPtr que = IntPtr.Zero;
            IntPtr prg = IntPtr.Zero;
            IntPtr kel = IntPtr.Zero;

            var pids = Introp.GetPlatformIds();
            //if (pids == null || pids.Length == 0)
            //{
            //    Console.WriteLine("zero platform");
            //    Console.ReadLine();
            //    return;
            //}
            var dic = new Dictionary<ulong, ulong[]>();
            for (int i = 0; i < pids.Length; i++)
            {
                dic.Add(pids[i], Introp.GetDeviceIds(pids[i]));
            }
            //foreach (var item in dic)
            //{
            //    Console.WriteLine(Introp.GetPlatformName(item.Key));
            //    if (item.Value == null || item.Value.Length == 0)
            //        continue;
            //    for (int i = 0; i < item.Value.Length; i++)
            //    {
            //        Console.WriteLine("\t" + Introp.GetDeviceName(item.Value[i]));
            //    }
            //}
            pid = dic.Keys.ToArray()[0];
            did = dic[pid][0];
            Console.WriteLine("选择的平台：" + Introp.GetPlatformName(pid));
            Console.WriteLine("选择的设备：" + Introp.GetDeviceName(did));
            ctx = Introp.CreateContext(did);
            que = Introp.CreateCommandQueue(did, ctx);
            string err;
            //使用cl文件
            //prg = Introp.CreateProgramFile(did, ctx, file, out err);
            //使用内置字符串
            prg = Introp.CreateProgramCode(did, ctx, CL_MulMatrix_Double, out err);
            kel = Introp.CreateKernal(prg, "matrix_mult");

            IntPtr pA = Introp.CreateClInputBuffer(ctx, A);
            IntPtr pB = Introp.CreateClInputBuffer(ctx, B);
            IntPtr pC = Introp.CreateClOutputBuffer(ctx, C);

            Introp.SetKernalParam(kel, 0, Ndim);
            Introp.SetKernalParam(kel, 1, Mdim);
            Introp.SetKernalParam(kel, 2, Pdim);
            Introp.SetKernalParam(kel, 3, pA);
            Introp.SetKernalParam(kel, 4, pB);
            Introp.SetKernalParam(kel, 5, pC);
            Int64[] glbal = new Int64[]
            {
                Ndim,Mdim
            };

            IntPtr pevent = Introp.Exec_CommandQueue(que, kel, 2, glbal);

            Introp.Exec_CommandQueueReadBuffer(que, pC, C);

            Introp.ReleaseMem(pA);
            Introp.ReleaseMem(pB);
            Introp.ReleaseMem(pC);
            Introp.ReleaseCommandQueue(que);
            Introp.ReleaseProgram(prg);
            Introp.ReleaseContext(ctx);

            return C;
        }

        #endregion GPU矩阵乘法

        #endregion 矩阵乘法

        #region 数学函数

        public const string CL_Cos_Float = @"
__kernel void cos_float(
	const int Ndim,
	__global const float* A,
	__global float* C)
{
	int i = get_global_id(0);
	//int j = get_global_id(1);
    float q = 0.0;

    C[i * Ndim + 0] = i;
    for(int t = 1;t<Ndim;t++)
    {
        switch(i)
        {
            case 0:
                C[i * Ndim + t] = sqrt(A[i * Ndim + t]);
            break;

            case 1:
                C[i * Ndim + t] = sin(A[i * Ndim + t])+100.0f;
            break;

            case 2:
                C[i * Ndim + t] = step(A[i * Ndim + t],-100.0f,-50.0f);
            break;

            case 3:
                C[i * Ndim + t] = tan(A[i * Ndim + t]);
            break;

            case 4:
                C[i * Ndim + t] = step(4 ,A[i * Ndim + t]);
            break;

            case 5:
                {
                    //C[i * Ndim + t] = mix(1000.0,10.0,0.3);
                    //C[i * Ndim + t] = sign(A[i * Ndim + t] - 4.0);
                    q = A[i * Ndim + t];
                    C[i * Ndim + t] = smoothstep(11,15,q);
                }
            break;

            default:
                C[i * Ndim + t] = A[i * Ndim + t];
            break;
        }
    }
}
";

        public static float[] Cos_Float(float[] A, int Ndim, int exec)
        {
            var C = new float[A.Length];
            Parallel.For(0, exec, i =>
            //for (int i = 0; i < exec; i++)
            {
                C[i * Ndim + 0] = i;
                for (int t = 1; t < Ndim; t++)
                {
                    switch (i)
                    {
                        case 0:
                            C[i * Ndim + t] = (float)Math.Sqrt(A[i * Ndim + t]);
                            break;

                        case 1:
                            C[i * Ndim + t] = (float)Math.Sin(A[i * Ndim + t]);
                            break;

                        case 2:
                            C[i * Ndim + t] = (float)Math.Cos(A[i * Ndim + t]);
                            break;

                        case 3:
                            C[i * Ndim + t] = (float)Math.Tan(A[i * Ndim + t]);
                            break;

                        case 4:
                            C[i * Ndim + t] = A[i * Ndim + t] < 4 ? 0 : 1;
                            break;

                        default:
                            C[i * Ndim + t] = A[i * Ndim + t];
                            break;
                    }
                }
            }
            );
            return C;
        }

        public static float[] Cos_Float_GPU(float[] A, int Ndim, int exec)
        {
            var C = new float[A.Length];

            ulong pid = 0;
            ulong did = 0;
            IntPtr ctx = IntPtr.Zero;
            IntPtr que = IntPtr.Zero;
            IntPtr prg = IntPtr.Zero;
            IntPtr kel = IntPtr.Zero;

            var pids = Introp.GetPlatformIds();
            var dic = new Dictionary<ulong, ulong[]>();
            for (int i = 0; i < pids.Length; i++)
            {
                dic.Add(pids[i], Introp.GetDeviceIds(pids[i]));
            }
            pid = dic.Keys.ToArray()[0];
            did = dic[pid][0];
            ctx = Introp.CreateContext(did);
            que = Introp.CreateCommandQueue(did, ctx);
            string err;
            //使用内置字符串
            prg = Introp.CreateProgramCode(did, ctx, CL_Cos_Float, out err);
            if (err != null)
            {
                Trace.WriteLine("OpenCL程序错误：" + err);
                Introp.ReleaseCommandQueue(que);
                Introp.ReleaseProgram(prg);
                Introp.ReleaseContext(ctx);
                return null;
            }
            kel = Introp.CreateKernal(prg, "cos_float");

            IntPtr pA = Introp.CreateClInputBuffer(ctx, A);
            //IntPtr pB = Introp.CreateClInputBuffer(ctx, B);
            IntPtr pC = Introp.CreateClOutputBuffer(ctx, C);

            Introp.SetKernalParam(kel, 0, Ndim);
            Introp.SetKernalParam(kel, 1, pA);
            Introp.SetKernalParam(kel, 2, pC);
            Int64[] glbal = new Int64[]
            {
                exec //,Mdim
            };

            IntPtr pevent = Introp.Exec_CommandQueue(que, kel, glbal);

            Introp.Exec_CommandQueueReadBuffer(que, pC, C);

            Introp.ReleaseMem(pA);
            //Introp.ReleaseMem(pB);
            Introp.ReleaseMem(pC);
            Introp.ReleaseCommandQueue(que);
            Introp.ReleaseProgram(prg);
            Introp.ReleaseContext(ctx);

            return C;
        }

        private static void TestMath()
        {
            Console.Title = "矩阵计算测试";
            Console.WriteLine("数学公式测试");

            int Ndim = 8;
            int exec = 8;
            int szA = Ndim * exec;
            float[] A = new float[szA];
            SpecMatrix(A, Ndim, exec);
            var c = Cos_Float(A, Ndim, exec);
            var gc = Cos_Float_GPU(A, Ndim, exec);
            PrintMatrix(A, Ndim, exec);
            Console.WriteLine("\t---CPU---");
            PrintMatrix(c, Ndim, exec);
            Console.WriteLine("\t---GPU---");
            PrintMatrix(gc, Ndim, exec);
        }

        #endregion 数学函数

        #region 同步机制

        public const string CL_Sync = @"
__kernel void cos_float(
	const int Ndim,
	__global const float* A,
	__global float* C)
{
	int i = get_global_id(0);
	//int j = get_global_id(1);

//内置的工作组内同步
work_group_barrier(CLK_LOCAL_MEM_FENCE) //
work_group_barrier(CLK_GLOBAL_MEM_FENCE)
work_group_barrier(CLK_IMAGE_MEM_FENCE)

    C[i * Ndim + 0] = i;
    for(int t = 1;t<Ndim;t++)
    {
        switch(i)
        {
            case 0:
                C[i * Ndim + t] = sqrt(A[i * Ndim + t]);
            break;

            case 1:
                C[i * Ndim + t] = sin(A[i * Ndim + t]);
            break;

            case 2:
                C[i * Ndim + t] = cos(A[i * Ndim + t]);
            break;

            case 3:
                C[i * Ndim + t] = tan(A[i * Ndim + t]);
            break;

            default:
                C[i * Ndim + t] = A[i * Ndim + t];
            break;
        }
    }
}
";

        #endregion 同步机制

        #region 全局内存

        private const string CL_Global_Memory_Array = @"
typedef struct AStructure{
  float a;
  float b;
}AStructure;
__kernel
void aFunction(__global AStructure *inputOutputBuffer)
{
    __global AStructure *inputLocation = inputOutputBuffer + get_global_id(0);
    __global AStructure *outputLoacation = inputOutputBuffer + get_global_size(0) + get_global_id(0);
    outputLoacation->a = inputLocation->a * -1;
    outputLoacation->b = (*inputLocation).b + 3.0f;
}
";

        #endregion 全局内存

        #region 局部内存

        private const string CL_LocalMem = @"
__kernel
void localAccess(__global float *A, __global float *B, __local float *C)
{
    __local float aLocalArray[1];

    if (get_local_id(0) == 0)
    {
        aLocalArray[0] = A[0];
    }
    C[get_local_id(0)] = A[get_global_id(0)];

    work_group_barrier(CLK_LOCAL_MEM_FENCE);

    float neighborSum = C[get_local_id(0)] + aLocalArray[0];
    if (get_local_id(0) > 0)
    {
        neighborSum = neighborSum + C[get_local_id(0) - 1];
    }
    B[get_global_id(0)] = neighborSum;
}";

        public static float[] CL_Local_Memory(float[] A, float[] B, int Ndim, int Mdim)
        {
            var ret = new float[B.Length];

            ulong pid = 0;
            ulong did = 0;
            IntPtr ctx = IntPtr.Zero;
            IntPtr que = IntPtr.Zero;
            IntPtr prg = IntPtr.Zero;
            IntPtr kel = IntPtr.Zero;

            var pids = Introp.GetPlatformIds();
            var dic = new Dictionary<ulong, ulong[]>();
            for (int i = 0; i < pids.Length; i++)
            {
                dic.Add(pids[i], Introp.GetDeviceIds(pids[i]));
            }
            pid = dic.Keys.ToArray()[0];
            did = dic[pid][0];
            ctx = Introp.CreateContext(did);
            que = Introp.CreateCommandQueue(did, ctx);
            string err;
            //使用内置字符串
            prg = Introp.CreateProgramCode(did, ctx, CL_LocalMem, out err);
            if (prg == IntPtr.Zero)
            {
                Console.WriteLine(err);
                return ret;
            }
            kel = Introp.CreateKernal(prg, "localAccess");

            IntPtr pA = Introp.CreateClInputBuffer(ctx, A);
            IntPtr pB = Introp.CreateClInputOutputBuffer(ctx, B);

            Introp.SetKernalParam(kel, 0, pA);
            Introp.SetKernalParam(kel, 1, pB);
            Introp.SetKernalLocalParam(kel, 2, Math.Max(Ndim, Mdim));

            Int64[] glbal = new Int64[]
            {
                Ndim*Mdim
            };

            IntPtr pevent = Introp.Exec_CommandQueue(que, kel, glbal);

            Introp.Exec_CommandQueueReadBuffer(que, pB, ret);

            Introp.ReleaseMem(pA);
            Introp.ReleaseMem(pB);
            Introp.ReleaseCommandQueue(que);
            Introp.ReleaseProgram(prg);
            Introp.ReleaseContext(ctx);

            return ret;
        }

        public static void TestLocalMemory()
        {

            Console.Title = "局部存储测试";
            Console.WriteLine("局部存储测试");

            int Ndim = 6;
            int Mdim = 5;
            int szA = Ndim * Mdim;
            float[] A = new float[szA];
            float[] B = new float[szA];
            NumberMatrix(A, Ndim, Mdim);
            var c = CL_Local_Memory(A, B, Ndim, Mdim);
            PrintMatrix(A, Ndim, Mdim);
            Console.WriteLine("\t---GPU---");
            PrintMatrix(c, Ndim, Mdim);
        }

        #endregion 局部内存

        private static void Main()
        {
            GPUEnvironment.Init();

            GPUEnvironment.Platforms.ToList().ForEach(p =>
            {
                Console.WriteLine(p.ToStringEx());
                p.Devices.ToList().ForEach(d =>
                {
                    Console.WriteLine(d.ToStringEx());
                });
            });

            TestMath();

            TestLocalMemory();

            TestMulMatrix();

            Console.WriteLine("...");
            Console.ReadLine();
        }
    }
}